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^sO We propose a new approach to finite density QCD based on a histogram method with phase 

quenched simulations at finite chemical potential. Integrating numerically the derivatives of 

{NJ the logarithm of the quark determinant with respect to the chemical potential, we calculate the 

i— ( reweighting factor and the complex phase of the quark determinant. The complex phase is han- 

dled with a cumulant expansion to avoid the sign problem. We examine the applicability of this 

t— I method. 
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1. Introduction 

In order to reveal the phase structure of QCD, which is relevant to the study of the the early 
universe, the core of the neutron star and the heavy ion collisions, it is indispensable to study QCD 
by first principle lattice simulations. The lattice simulations, however, have the notorious sign prob- 
lem at non-zero quark chemical potential }i. Although there have been proposed several approaches 
to study finite density QCD, a satisfactory method reliable at large quark chemical potential is still 
lacking. We propose a new approach to finite density QCD by means of the histogram method and 
the reweighting technique [1,2] together with phase quenched simulations, in which the Monte 
Carlo ensemble is generated without the complex phase of the quark determinant. The complex 
phase is handled with a cumulant expansion to evade the sign problem. In this report, we examine 
the applicability of our method; in particular, we discuss the overlap problem and the convergence 
of the cumulant expansion. 

2. Histogram method 

In order to calculate thermodynamic quantities such as the pressure, we need to calculate the 
expectation values of the plaquette and the quark determinant. If we are interested in quantities 
which depend only on them, the histogram method enables us to evaluate the expectation values 
from the probability distribution function of the plaquette and the quark determinant. Here, we 
discuss the case of the degenerate Nf flavor case. An extension to the non-degenerate case is 
straightforward. We label the gauge configurations by the space-time averaged plaquette P and 
the absolute value of the quark determinant, F{[i) = Nf\n |detM(/i)/ detM(0)|. Decomposing the 
quark determinant as (detM( J u)) A ' f = e' 6 ^ \detM(n)\ Nl , the partition function normalized at zero 
chemical potential can be written as 

W$ = W0j I ® Uei9W \ d ^)\ Nf ^ P = J dPdF (e idW ) {pF Wo(P,F,l5,^ 

(2.1) 

where the probability distribution function, 

w,{P\F\^^) = J ^J&U8{P'-P[U})8{F'-F[U})\dttM{ l i)\^ (2.2) 

is obtained by the histogram of P and F in the phase quenched simulation, and 

m \ _ J We ie > 8 {P' - P[U] ) 8 {F' -F[U])\ detM(n)\ Nf e 6 ^ p 

/(p',f>) " f@U8(P'-P[U])8{Fi-F[U})\detM(n)\ Nf e 6 P N ^ p 
((e*MS(P> -P[U])S(F> -F[U]))) (M 

((8(P'-P[U])8(F'-F[U]))) {M (23) 

is the expectation value of the complex phase of the quark determinant with fixed P' and F'. The 
double bracket indicates the expectation value in the phase quenched simulation. Here, ./V s ; te = 
A^ 3 x N t is the number of lattice sites and /3 =6/g 2 . Note that (e' d ) does not depend on j8 since we 
can factor out e 6 ^ NsiteP from both the numerator and the denominator in Eq. (2.3). Introducing the 
effective potential Vb = — lnwo, the ratio of the partition function, Eq. (2.1), can be written as 

= f dPdFe-i^M-^^} = f dPdFe-WM, V = V -ln(e w ^). (2.4) 
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In order to calculate the partition function precisely, 
we have to evaluate the P and F integration in Eq. (2.4) 
accurately: namely, we have to generate gauge config- 
urations near the minimum of V, which dominates the 
integral in Eq. (2.4). The histogram of P and F at a sin- 
gle simulation point (j6,/i) provides the effective poten- 
tial Vo covering a limited region in the (P,F) plane. If 
the complex phase (2.3) has a large P and/or F depen- 
dence, the minimum of V differs from that of Vo, and we 
do not have sufficient number of configurations near the 
minimum of V (see Fig. 1). In such a case, there is no 
sufficient overlap between the important region for the 

integral of Eq. (2.4) and the region where Vo is evaluated precisely by measuring the histogram in 
the phase quenched simulation. 

This overlap problem can be circumvented by combining the histograms at several simulation 
points with the aid of the reweighting method. The effective potential Vo at (j8,jti) can be obtained 
by that at (J3o,jUo) as follows: 



>P(F) 

Figure 1: A schematic figure for the 
overlap problem. 



V Q (P,F,p,n) = V (P,F, ft, Mo) - lnR(P,F, p , ft, fi, juo). 
Here /?(P',F',j8,j8o,jU,jUo) = w (P,F, j3 , /w (P,F, ft, Mo) is the reweighting factor, 



(2.5) 



R(P , ,F',p,p ,li,^)=e 



6(P-p )N sltB P' 



8(P'-P[U])8(F'-F[U}) 



detM(p) 
detMQxo) 



(ft,AK>) 



((S(P>-P[U])S(F>-F[U]))) 



(ft, Mo) 



(2.6) 



Determination of R is reduced to that of the expectation value of the quark determinant in the 
phase quenched simulation when /3 = ft. Under a /3 shift with keeping = Mo, the slope of the 
effective potential changes by a constant factor while the curvature remains the same. Combining 
the effective potentials Vo at various simulation points using Eq. (2.5), we can have an enough 
overlap with the minimum of V even if the complex phase has a large P and/or F dependence. We 
can thus avoid the overlap problem. 

Because the phase quenched simulations in two-flavor QCD correspond to the case of isotriplet 
chemical potentials, a comment is in order about the influence of the pion condensed phase. The 
large isotriplet chemical potential induces the pion condensation [3]. In the pion condensed phase 
( e ' 9 )(p,F) * s ex P ecte d to vanish as has been suggested in model calculations [4, 5]. This implies 
that Vo(P,F) and V(P,F) = Vo(P,F) — ln(e~' e ^ pF ^ have no overlap inside the condensed phase 
and the partition function Z(/3,m) is dominated by configurations outside the condensed phase. 
Thus, we do not need to generate configurations with (e' e )( PF j = 0, which have no contribution to 
the integral in Eq. (2.1). 



3. Cumulant expansion for the complex phase of the quark determinant 

Even though we can circumvent the overlap problem, a large fluctuation of the phase of the 
quark determinant at large chemical potential leads to a frequent change of the sign of the complex 
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phase. In this case, Monte Carlo simulations suffer from the sign problem. We exploit the cumulant 
expansion for the expectation value of the complex phase to avoid the sign problem, 



(e 



(p,F) = exp 



6- 



i 

3! 



'0- 



(3.1) 



The cumulants agree with the central moments up to the third order, while they differ at higher 
orders. For instance, the fourth order cumulant is given by 



e-(fl) 



(P-F) 



o -(e) 



(P,F) 



(3.2) 



(RF) \ 'I (P,F) 

The odd-order cumulants change the sign under the flip of the sign of the chemical potential, 
H -R- —pL, which transforms quarks into antiquarks. Accordingly, only the even-order cumulants 
survive if the system is invariant under this transformation, and the complex phase now becomes 



(e ie ^) {P , F) =^V 



(3.3) 



We stress that the right-hand side is real and positive once we drop the odd-order cumulants from 
the symmetry. By applying the cumulant expansion for the complex phase, the sign problem is 
reduced to the convergence problem of the cumulant expansion; namely, we have no sign problem 
if the cumulant expansion converges. 

An ideal case for the convergence of the cumulant expansion is the case that the phase 6 has 
a Gaussian distribution. In such a case, only the second-order cumulant survives, (e' 9 ^)( P ^ = 

exp — (0 2 )( PF ) /2 - If we calculate 0(jU) =./Vf3m[lndetM( i u)] in the limited range [—71, 7t) taking 
into account the periodicity of the complex phase (e' e ), the phase distribution may have no resem- 
blance to the Gaussian distribution. It is essential for the convergence of the cumulant expansion 
to calculate the phase of the quark determinant such that the distribution takes of nearly a Gaussian 
form. In this study, we define the phase in the range — oo < < oo. Instead of calculating detM(/i) 
directly, we measure the the derivatives of lndetM(/x) with respect to ju, and then calculate the 
phase of the quark determinant by integrating the derivatives over ju, 



e(jU) =Af f 3m[lndetM( i u)] =N f / 3m 

Jo 



5(lndetM(^i)) 



(3.4) 



Note that this is not a Taylor expansion; namely, it is applicable to any values of the chemical 
potential. Conventional phase in the range [—71, 7i) is recovered by taking the principal value of 6 
with the period of 2ti. The integrand in Eq. (3.4), 



d(lndetMGu)) 
d(ji/T) 



Tr M _1 (/x) 



d(n/T)J 



(3.5) 



can be regarded as the sum of the local density operator denned at each lattice site. If the density 
operator has a small correlation length compared to the spatial size of the system, we expect the 
Gaussian distribution for the operator due to the central limit theorem [1]. The volume dependence 
of the convergence of cumulant expansion has been discussed in [6]. 

The absolute value of the quark determinant, which is used to label the gauge configurations, 
and the ratio of the quark determinant, which is needed to evaluate the reweighting factor Eq. (2.6), 
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can be also obtained by integrating the real part of the derivatives without further computational 

costs, 



F(ji) = Nfln 
C(ji) = Nfln 



detM(n) 



detAf(O) 
detM(ju) 



detM(jUo) 



Jo 
Jiki 



d(lndetM(/i))" 
WJT) _ 
d(lndetM(») 

wm 



(3.6) 
(3.7) 



We note that 6(fl), F(fi), and C(ju) can be obtained as continuous functions of pL in this approach. 
In addition, the statistical errors for the reweighting factor R are expected to be small for fixed F 
since F and C are strongly correlated. 

4. Numerical simulations and the results 

In this study, we use the RG-improved Iwasaki action for gauge action and the Nf = 2 0(a)- 
improved Wilson quark action with csw = (1 — 0.8412/J -1 ) -3 / 4 . The ratio of pseudoscalar and 
vector meson masses at T = }X = are set to raps /ray = 0.8. We generate gauge configurations on 
8 3 x 4 lattice with the complex phase of the quark determinant removed. The measurement of the 
first and the second derivatives of IndetM (ju) with respect to fJL, which are used to interpolate the 
integrands in Eq. (3.4), (3.6) and (3.7), has been done every 10 trajectories. We employ the random 
noise method of [6] with 50 noises. The statistics is the order of 0(1000). 
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Figure 2: The effective potential Vo(F) at \ljT = 2.0 evaluated at three different simulation points. 

The effective potential without the complex phase, Vo(F) = — \awo(F), at \ijT = 2.0 evalu- 
ated at three different simulation points, (fi,Ho/T) = (1.5, 1.6), (1.5,2.0), (1.5,2.4), is drawn in 
Fig. 2. Vo(F) is normalized such that Vq(F) = at the minimum for each simulation. We observe 
that the three data sets covering different ranges nicely fall on one curve. Although Vn at the single 
simulation point (fi,Ho/T) = (1.5,2.0) covers only the narrow range centered around F ~ 50, the 
effective potential obtained from the histograms at different simulation points by the use of the 
reweighting method widely range in the F direction. Moreover, we see that the statistical errors of 
the effective potential, which stem from the reweighting factor, are small as we expected. 

The distribution of the phase of the quark determinant is depicted in Fig. 3. The dashed curves 
are the fitted results with a Gaussian function. Bf is the fourth-order Binder cumulant normalized 
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Figure 3: The distribution of the phase of the quark determinant at ji/T = 0.4 (upper panels) and ji/T = 2.4 
(lower panels). The dashed curves are the fitted results with the Gaussian function. B® is the fourth-order 
Binder cumulant normalized such that B® = 3 for the Gaussian function. 



suchthatB^ = 3 for a Gaussian distribution, i.e.,^ = (6 4 ) c /(d 2 )l + 3. The upper and lower panels 
are the results at /J./T = 0.4 and 2.4, respectively. We see that the phase distribution gets broader 
as the chemical potential increases. Furthermore, at large chemical potential (lower panels), we 
observe that the phase distributes broadly at low temperature. The important point is that the phase 
distribution evaluated by Eq. (3.4) can be well approximated by a Gaussian function even in the 
high density region fx/T > 1. This promises a good convergence in the cumulant expansion for the 
complex phase of the quark determinant. 

Fig. 4 shows the second and the fourth order of the cumulants as a function of F (upper panels) 
and P (lower panels) at fl/T = jx /T = 0.4 (left panels) and n/T = Ho/T = 1 .2 (right panels). We 
observe that the second order cumulant increases with \ijT both as a function of F and P. On the 
other hand, the fourth order cumulant is consistent with zero within the statistical errors although 
the eiTors grow with fJ./T. We do not see a clear F and P dependence of the second and fourth 
cumulants in these parameter region. 

5. Summary 

In this study, we proposed a new approach to finite density QCD based on the histogram 
method and the reweighting technique with phase quenched simulations. We apply the cumulant 
expansion for the complex phase of the quark determinant. We found that the reweighting technique 
combined with the histogram method works well and we obtained the effective potential covering a 
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Figure 4: The second and the fourth order of the cumulants as a function of F (upper panels) and P (lower 
panels) at \ljT = Ho/T = 0.4 (left panels) and n/T = /Xq/T = 1.2 (right panels). 



wide range. Moreover, we proposed a way to calculate the phase 6 of the quark determinant, which 
is not constrained in the range [—71, 7l). We showed that the distribution of 6 becomes wide as Li/T 
increases and the phase distribution can be well approximated by a Gaussian function both at small 
and large chemical potential. The second order cumulant increases with the chemical potential, 
while the fourth order cumulant is consistent with zero although the statistical error increases with 
jU. A comprehensive analysis in a wide parameter region is in progress. 

This work is in part supported by Grants-in-Aid of the Japanese Ministry of Education, Cul- 
ture, Sports, Science and Technology, (Nos.20340047, 21340049, 22740168, 23540295) and by the 
Grant-in-Aid for Scientific Research on Innovative Areas (Nos. 20105001, 20105003, 2310576). 
HO is supported by the Japan Society for the Promotion of Science for Young Scientists. 

References 

[1] S. Ejiri, Phys. Rev. D 11 (2008) 014508, arXiv:0706.3549. 

[2] H. Saito et al. (WHOT-QCD Collaboration), Phys. Rev. D 84 (2011) 054502, arXiv: 1 106.0974. 
[3] D. Son and M. A. Stephanov, Phys. Rev. Lett. 86 (2001) 592, hep-ph/0005225. 
[4] J. Han and M. Stephanov, Phys. Rev. D 78 (2008) 054507, arXiv:0805.1939. 
[5] Y. Sakai et al, Phys. Rev. D 82 (2010) 096007, arXiv: 1005.0993. 

[6] S. Ejiri et al. (WHOT-QCD Collaboration), Phys. Rev. D 82 (2010) 014508, arXiv:0909.2121. 



7 



